##### Replicates the analysis in the frontiers article.  Note that
##### bootstrapped CIs will be slightly different on every run.  All
##### page numbers are for word doc and may not match typeset version.

## Load mosaic library
library(mosaic)

#### Experiment 1

## Prep
rm(list=ls())
height <- read.csv("frontiers1.csv")

## Calculate standardized deviations from true value
height$ErrorRatio <- with(height, (TrueHeight-Estimate)/TrueHeight)

### T-test on page 5 (result for exp 1)

## Helper function, to make boots easier
my.t.test <- function (dat) {
  tst <- t.test(ErrorRatio ~ Condition, data=dat)
  var.down = with(dat, var(ErrorRatio[Condition=="down"]))
  var.up = with(dat, var(ErrorRatio[Condition=="up"]))
  n.down = sum(dat$Condition=="down")
  n.up = sum(dat$Condition=="up")
  sd.all = sqrt(((n.down-1)*var.down + (n.up-1)*var.up)/(n.down+n.up-2))         
  d = (tst$estimate[1] - tst$estimate[2]) / sd.all
  return(list(tst=tst, d=d))
}

tst.res <- my.t.test(height)
tst.res
tst.boot <- do(1000)*my.t.test(resample(height))
quantile(unlist(tst.boot$d), probs=c(0.025, .05, .95, .975))

#### Experiment 2

## Clear workspace, prep
rm(list=ls())
bos <- read.csv("frontiers2.csv")
bos <- subset(bos, Drop==0) # pair both given same sheet
bos$Room[bos$Room==1] <- "up"
bos$Room[bos$Room==2] <- "down"
bos$Choice <- as.factor(bos$Choice)
bos$Condition <- as.factor(bos$Condition)

### pg 7 test of pattern choice different from 0.5
tst <- prop.test(bos$Pattern=="A")
tst

## Cohen's w (equals phi here) for above test, w/bootstrapped CIs
sqrt(tst$statistic/nrow(bos))
test.boot <- do(1000)*prop.test(subset(resample(bos), Condition!=3)$Pattern=="A")$statistic
cohen.w.boot <- unlist(sqrt(test.boot/nrow(bos)))
quantile(cohen.w.boot[1:10], probs=c(0.025, .05, .95, .975))

### Table 2
fi <- glm(Choice=="3" ~ Room*Condition, data=bos, family="binomial")
summary(fi)

### Pairwise analysis

## Set up pairs.  Assumes pairs in consecutive rows in data.
pairs <- bos[seq(1, nrow(bos), 2),]  # take first person in each pair
pairs <- subset(pairs, Winnings > 0) # drop non-coordination cases
## Now code as TRUE if down wins, FALSE if up wins
pairs$Downwin <- with(pairs, (Room=="down" & Winnings=="3") | (Room=="up" & Winnings=="2"))

## Test of proportion equality convenience function
between.test <- function (exp1, exp2) {
  ## One-tailed test, no small sample correction
  prop.test(x=c(sum(exp1$Downwin), sum(exp2$Downwin)),
            n=c(nrow(exp1), nrow(exp2)), 
            alternative="g", correct=FALSE)
}

### Pairwise test, supporting logistic regression at individual level
btwn<-between.test(subset(pairs, Condition=="1"),  subset(pairs, Condition=="2"))
btwn

## Cohen's h for above test, w/bootstrapped CIs
2*(asin(sqrt(btwn$estimate[1]))-asin(sqrt(btwn$estimate[2])))
btwn.boot <- do(1000)*between.test(resample(subset(pairs, Condition=="1")),  resample(subset(pairs, Condition=="2")))$estimate
cohen.h.boot <- 2*(asin(sqrt(btwn.boot[,1]))-asin(sqrt(btwn.boot[,2])))
quantile(cohen.h.boot, probs=c(0.025, .05, .95, .975))